######################################################################################################
# ADICA - Automatic Differentiation Integrated Catch Analysis
#
# Version 0.12 26/02/2009 10:57:21
#
# Author: Mark Payne, mpa@aqua.dtu.dk
# DTU-Aqua, Charlottenlund, DK
# 
# An AD Model Builder Implementation of ICA, with an interface to FLR
#
# Developed with:
#   - R version 2.8.0
#
# Changes:
#
# To be done:
#   *  Settle on a format for values to exclude. NA or -1? Or do the exclusion actively
#      in the R code?
#   *  Index model slot and index type should be enumerated
#
# Checks on input data:
#   *  Each FLIndex has a correctly set type attribute
#   *  Require at least four ages. Or more? One of which is a plus group
#   *  Catch slot is occupied
#   *  Stock plusgroup is set
#   *  Control object sizes must be consistent with the stock and indices object. Cross check if needed
#   *  Index ages must be continuous
#
# This script is subject to Mark's version of the "BEER-WARE" LICENSE.
# Mark Payne wrote this file. It is not intended for use by other authors.
# If it breaks, it is not my fault. Support is not provided. But if you are determined
# to use it anyway, and we meet some day, then, if you think this stuff is worth it,
# you can buy me a beer in return. Mark.
####################################################################################################

### ======================================================================================================
### Initialise system
### ======================================================================================================
# Start with house cleaning
rm(list = ls(all.names=TRUE)); gc(); graphics.off()
ver <- "\nADICA v 0.12\n"; cat(ver)
ver.datetime   <- "26/02/2009 10:57:21\n\n";
cat(ver.datetime); start.time <- proc.time()[3]
options(stringsAsFactors=FALSE)
library(FLICA)


write.adica.dat <- function(stck,tun,ctrl,output.file) {
### ======================================================================================================
### Setup data file
### ======================================================================================================
old.opt <- options("scipen"=100)
cat("#######################################\n# Automatic Differentiation ICA (ADICA)",
    "\n# This file is automatically generated","\n# by the corresponding R script.",
    "\n#######################################",file=output.file)

### ======================================================================================================
### Write Stock info to dat file
### ======================================================================================================
cat("\n#\n##########################\n# Stock Information\n##########################\n",file=output.file,append=TRUE)
#Write the components of the file
cat("# stckRange: Stock range information: minage maxage plusgroup minyear maxyear startf endf\n",file=output.file,append=TRUE)
write.table(stck@range,file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)
cat("# nAges: Number of ages\n",dims(stck)$age,"\n",file=output.file,append=TRUE)
cat("# nYears: Number of years\n",dims(stck)$year,"\n",file=output.file,append=TRUE)
cat("# canum: Catch at age in numbers\n",file=output.file,append=TRUE)
write.table(stck@catch.n@.Data,file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)
cat("# west: Weight in the stock\n",file=output.file,append=TRUE)
write.table(stck@stock.wt@.Data,file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)
cat("# natM: Natural Mortality\n",file=output.file,append=TRUE)
write.table(stck@m@.Data,file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)
cat("# propMat: Proportion Mature\n",file=output.file,append=TRUE)
write.table(stck@mat@.Data,file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)
cat("# fProp: Proportion of F occuring before spawning\n",file=output.file,append=TRUE)
write.table(stck@harvest.spwn@.Data,file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)
cat("# mProp: Proportion of natural Mortality occuring before spawning\n",file=output.file,append=TRUE)
write.table(stck@m.spwn@.Data,file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)

### ======================================================================================================
### Write Indices to dat file
### ======================================================================================================
#Get number of indices
noIdxs  <- length(tun)
#Extract data into a list format. Set NAs to negative numbers.
idxDat <- list()
for(i in 1:noIdxs) {
                idx <- tun[[i]]
                dat <- data.frame(Index=i,expand.grid(dimnames(idx@index)[1:2]),index=as.vector(idx@index@.Data),index.var=as.vector(idx@index.var@.Data))
                idxDat[[i]] <- split(dat,list(dat$age))
            }
idxDat    <- do.call(c,idxDat)
idxDat    <- lapply(idxDat,function(x) {
                if(unique(x$age)=="all") {x$age <- -1}
                return(x)})
nIDs      <- length(idxDat)
idxLength <- sapply(idxDat,nrow)
#Extract range and index type data
idxsRange  <- t(sapply(tun,range))
idxsType   <- factor(sapply(tun,type),levels=c("number","ssb"))
#Now write the indices component of the file
cat("#\n#\n##########################\n# Indices\n##########################\n",file=output.file,append=TRUE)
cat("# noIdxs : Number of indices\n",noIdxs,"\n",file=output.file,append=TRUE)
cat("# nIDs : Number of unique survey/age combinations\n",nIDs,"\n",file=output.file,append=TRUE)
cat("# idxsRange: Index range information: min max plusgroup minyear maxyear startf endf\n",file=output.file,append=TRUE)
write.table(idxsRange,file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)
cat("# idxTyp : Type of each index (1=number, 2=ssb) \n",idxsType,"\n",file=output.file,append=TRUE)
cat("# nIdxPts : Total number of index data items (rows) \n",sum(idxLength),"\n",file=output.file,append=TRUE)
for(i in 1:nIDs){
    cat("# idxIn (ID, indexNumber, age, year, index, index.var)\n",file=output.file,append=TRUE)
    write.table(cbind(i,idxDat[[i]]),file=output.file,na="-1",append=TRUE,row.names=FALSE, col.names=FALSE, quote=FALSE)
}

### ======================================================================================================
### Write Control info to dat file
### ======================================================================================================
cat("#\n##########################\n# ICA Control parameters",
    "\n# index.model types (0=absolute,1=linear,2=power)",
    "\n# sr (0=FALSE,1=TRUE)\n##########################",file=output.file,append=TRUE)
lapply(slotNames(ctrl),function(slt) {
    cat("\n#",slt,"\n",as.numeric(slot(ctrl,slt)),file=output.file,append=TRUE)
})

}


write.adica.pin <- function(ica,pin.file) {
### ======================================================================================================
### Now write the results of the assessment as parameters
### ======================================================================================================
cat("# F_y:\n",file=pin.file)
write.table(t(subset(ica@param,Param=="F",select=6)),file=pin.file,quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)
cat("# Sel:\n",file=pin.file,append=TRUE)
write.table(t(subset(ica@param,Param=="Sel",select=6)),file=pin.file,quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)
cat("# TY_pop:\n",file=pin.file,append=TRUE)
write.table(t(subset(ica@param,Param=="TermN",select=6)),file=pin.file,quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)
cat("# LTA_pop:\n",file=pin.file,append=TRUE)
write.table(t(subset(ica@param,Param=="OldN",select=6)),file=pin.file,quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)
cat("# Recruitment Prediction:\n",file=pin.file,append=TRUE)
write.table(t(subset(ica@param,Param=="Rec",select=6)),file=pin.file,quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)
cat("# Qs:\n",file=pin.file,append=TRUE)
write.table(t(subset(ica@param,Param=="Q",select=6)),file=pin.file,quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)
cat("# Ks:\n",file=pin.file,append=TRUE)
write.table(t(subset(ica@param,Param=="K",select=6)),file=pin.file,quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)
cat("# SRR a and b:\n",file=pin.file,append=TRUE)
write.table(t(subset(ica@param,Param=="SRRa"|Param=="SRRb",select=6)),file=pin.file,quote=FALSE,row.names=FALSE,col.names=FALSE,append=TRUE)
}


### ======================================================================================================
### Write files
### ======================================================================================================
#herIIIa object
data(herIIIa)
source("FLICA Class.r")
ctrl <- FLICA.control(sep.nyr=5,sep.age=4,sep.sel=1.0,lambda.yr=rep(1,5),lambda.age=c(0.1,1,1,1,1,1,1,1,1),
                              lambda.sr=0,sr=FALSE,sr.age=0,index.cor=1)
ctrl@index.model <- ac(c(1,1,1))    #Use 0 for abs, 1 for linear, 2 for power - in that way, the value tells how many parameters to fit

write.adica.dat(herIIIa,herIIIa.tun,ctrl,"adica.dat")
write.adica.pin(herIIIa.ica,"adica.pin")

#North Sea herring object
load("NSH Assessment.RData")
NSH.tun[[1]]@type <- "ssb"
NSH.ctrl@index.model <- ac(c(2,1,1,1))
write.adica.dat(NSH,NSH.tun,NSH.ctrl,"adica.dat")
#write.adica.pin(NSH.ica,"adica.pin")
